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ABSTRACT 

We study numerically the evolution of an adiabatic relativistic fireball 
expanding into a cold uniform medium. We follow the stages of initial free 
expansion and acceleration, coasting and then deceleration and slowing down 
to a non-relativistic velocity. We compare the numerical results with simplified 
analytical estimates. We show that the relativistic self similar Blandford-McKee 
solution describes well the relativistic deceleration epoch. It is an excellent 
approximation throughout the relativistic deceleration stage, down to 7 ~ 5, 
and a reasonable approximation even down to 7 ~ 2 though the solution is 
rigorous only for 7 ^> 1. We examine the transition into the Blandford-McKee 
solution, and the transition from the solution to the non-relativistic self similar 
Sedov- Taylor solution. These simulations demonstrate the attractive nature of 
the Blandford-McKee solution and its stability to radial perturbations. 

Subject headings: hydrodynamics — relativity — shock waves — gamma-ray 
bursts 
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1. Introduction 

Sedov (1946), Taylor (1950) and Von Neumann (1947) discovered, in the forties, a 
self similar solution of the strong explosion problem, in which a large amount of energy 
is released on a short time scale in a small volume. This solution is known today as the 
"Sedov- Taylor" self similar solution. It describes a shock wave propagating into a uniform 
density surrounding. The shock wave and the matter behind it decelerate as more and more 
mass is collected. This solution describes well the adiabatic stage of a supernova remnant 
evolution. 

Blandford and McKee (1977) have later established a self similar solution describing 
the extreme relativistic version of the strong explosion problem. In this solution the Lorentz 
factor of the shock and the fluid behind it is much larger than unity. Such high Lorentz 
factors arise if the rest mass contained within the region where the energy E was released is 
much smaller than E. In other words when the region containing the energy is "radiation 
dominated" rather than matter dominated. Such a region was later termed a "fireball" . 

Cavallo and Rees (1978) have considered the physical processes relevant in a radiation 
dominated fireball as a model for gamma-ray bursts (GRBs). Goodman (1986) and 
Paczyhski (1986) have then considered the evolution of such a fireball. They have shown 
that a initially the radiation-pair plasma in a purely radiative fireball behaves like a fluid 
and it expands and accelerated under its own pressure until the local temperature drops 
to ~ 20keV, when the last pairs annihilate and the fireball becomes optically thin. Later 
Shemi and Piran (1990) considered matter contaminated fireballs. They have shown that, 
quite generally, all the initial energy will be transfered to the baryons in such fireballs whose 
final outcome is a shell of relativistic freely expanding baryons. Piran, Shemi and Narayan 
(1993) and Mesazaros Laguna and Rees (1993) have later carried these calculations in 
greater details. These works show that an initially homogeneous fireball will first accelerate 



-4- 



while expanding and then coast freely as all its internal energy was transformed to kinetic 
energy 

The surrounding matter will eventually influence the fireball after enough external 
matter has been collected and most of the energy has been transfered from the shell to the 
ISM. If the surrounding matter is diluted enough, then this will take place only after the 
initial free acceleration phase. This influence was considered by Mesazaros and Rees (1992), 
and Katz (1993), who suggested that the GRB is produced during this stage. The detailed 
shock evolution was later studies by Sari and Piran (1995). It is only after these stages, 
that the fireball have given the ISM most of its energy and then the self similar deceleration 
solution of Blandford-McKee applies. When the shock decelerates enough so that it is no 
longer relativistic, it is described by the Sedov- Taylor solution. 

Today it is widely accepted that GRBs involve relativistic expanding matter of this 
kind. While the GRB itself is produced, most likely, via internal shocks (Narayan, Paczyhski 
& Piran, 1992; Rees & Mesazaros, 1994; Sari & Piran, 1997). The observed GRB afterglow 
corresponds, on the other hand to the slowing down of this relativistic flow. This has 
led to an increasing interest in the fireball solution and in its various regimes. In this 
paper we study the evolution of a homogeneous fireball focusing on its interaction with the 
surrounding matter. We do not consider here internal shocks, which arise due to interaction 
within the relativistic flow and require nonuniform velocity. 

We have developed a spherically symmetric relativistic Lagrangian code based on a 
second order Gudnov method with an exact Riemann solver to solve the ultra-relativistic 
hydrodynamics problem. With this code, it is possible to track the full hydrodynamical 
evolution of the fireball within a single computation, from its initial "radiation dominated" 
stage at rest through its acceleration, coasting, and shock formation, up to the relativistic 
deceleration and finally to the Newtonian deceleration. With typical parameters this 
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computation spans more than eight orders of magnitudes in the size of the fireball and more 
than twenty orders of magnitudes in its density. We describe these computations here. We 
show that, quite generically, the solution converges during the relativistic deceleration phase 
to the Blandford-McKee solution and then it transforms to the Sedov- Taylor solution. Even 
though the attractive nature of the Blandford-McKee solution suggests that it is stable 
we explore explicitly the stability of this solution and we show that it is stable to radial 
perturbations. 

We review, first, in section |2] the current analytic understanding of the fireball evolution 
thorough the following stages: (i) free acceleration and coasting, (ii) energy transfer (iii) 
relativistic self similar solution and (iv) Newtonian Sedov- Taylor solution. We discuss the 
numerical results in section |3|. In section |] we examine the evolution of perturbations to the 
Blandford-McKee solution. We discuss the implications of these results in section |5|. 



2. Analytic Estimates 

The evolution of a fireball is characterized by several phases. The transitions between 
these phases are determined by several critical radii that are summarized in table |T]. 

Table 1: Critical Radii 



Ro 


Initial Fireball Size 




Rl 


Coasting 


RoV 


R s 


Spreading (and Internal Shocks) 


Rov 2 


Ra 


External Shocks (the NRS Case) 


/3 /4 A l/4 


R^y 


External Shocks (the RRS Case) 




Rn 


Relativistic Reverse Shock 




I 


Sedov Length 


(£/Vi) 1/3 
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2.1. Free Acceleration and Coasting 

We consider a homogeneous fireball of energy E and a baryonic load of total mass 
M confined initially in a sphere of radius Rq. We define the dimensionless entropy (or 
the initial random Lorentz factor) r) = E/M . This fireball expands into a surrounding 
low density medium (with a density pi) which we will refer to as the ISM. This can be 
considered to be a free expansion in its initial stage. After a short acceleration phase, 
the motion becomes highly relativistic. Conservations of baryon number, energy and 
momentum yield the following conservation laws along a null flow line of each fluid element 
in the shell (Piran, Shemi & Narayan, 1993): 

r 2 p7, r 2 p 3,/4 7 and r 2 (Ap + p)^y 2 = constant, (1) 

where r(t), j(t), p(t) and p(t) are the radius, Lorentz factor, pressure and rest mass 
density of the fluid element, respectively. In this paper, distance, time, velocity and 
the corresponding Lorentz factors are measured in the observer frame. Thermodynamic 
quantities (p and p) are measured in the local fluid frame. We use units in which the speed 
of light c = 1. The above equations assume an adiabatic gas index of 4/3. 

Initially the fireball is extremely hot (p ^> p), so that equation ([I]) yields 

7 oc r, p oc r~ 3 and p oc r -4 , (2) 

(Goodman, 1986; Paczyhski, 1986; Shemi & Piran, 1990). The fireball is approximately 
homogeneous in the local frame, but due to relativistic effects it appears it appears to an 
observer at rest as a narrow shell with a radial width A ~ r/7 ~ Rq (Shemi & Piran, 
1990; Piran, Shemi & Narayan, 1993). As the fireball expands, the internal energy is 
converted into kinetic energy of the baryons. At the radius Rl = r)R , the fireball uses up 
all the internal energy and the approximation p ^> p breaks down. This is the end of the 
acceleration phase. 
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Now, the internal energy of the fireball becomes negligible compared to the rest mass 
energy (p <C p), and equation (|l|) yields 

7 = constant, p oc r~ 2 and p oc r~ 8//3 , (3) 

(Piran Shemi & Narayan, 1993). The fireball behaves like a pulse of energy with a frozen 
radial profile propagating at almost the speed of light. 

2.2. Spreading, The Reverse Shock and Energy Transfer 

This frozen pulse approximation on which equation (P is based breaks down ultimately 
at the radius R s = Rq1] 2 . Each fluid shell moves with a slightly different velocity and the 
fireball begins to spread at R s . Internal shocks will take place around R s if the fireball is 
inhomogeneous and the velocity is not a monotonic function of the radius. As mentioned 
earlier, these shocks produce, most likely, the observed GRB. However, even under optimal 
condition they cannot convert more than about a quarter of the kinetic energy to radiation. 
Hence even an inhomogeneous shell will continue carrying ample kinetic energy beyond this 
stage. We consider here only homogeneous fireball. The possible effect of internal shocks on 
the fireball evolution has been discussed extensively in another papers (Kobayashi, Piran & 
Sari, 1997; Daigne & Mochkovitch, 1997) 

The coasting can also end if the surrounding matter begins to influence the shell. The 
interaction between the shell and the ISM can be described by two shocks: a forward shock 
propagating into the ISM and a reverse shock propagating into the shell. Sari and Piran 
(1995) have defined three critical radii in this respect: = / 3//2 /A 1//2 ^ 2 where the energy 
density produced by the shocks becomes high enough so that the reverse shock is relativistic 
and begins to reduce the Lorentz factor of the shell considerably; Ra = / 3//4 A 1/ ' 4 where the 
reverse shock crosses the shell; and R 1 = l/r) 2 ' 3 where the mass of the shocked ISM is Mq/t]. 
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Here / = (E / pi) 1 ^ 3 is the Sedov length. Fortunately, a simple relations between these four 
radii can be given in terms of the dimensionless variable £ = (//A) 1 / 2 ?? -4 / 3 ; 

eR s = v&a = Ri = Rn/Z- (4) 

If initially £ > 1 then R s is the smallest radius. The shell begins to spread at R s . After 
that the width A satisfies A = r/7 2 oc r and the scaling of the shell parameters becomes 

7 = constant, p oc r~ 3 and p oc r -4 . (5) 

During the spreading phase the value of £ > 1 decreases. However, as long as £ > 1 the 
relation R s < Ra < R y < Rn is valid. When £ ~ 1 these different radii become comparable: 
the reverse shock crosses the shell, it becomes mildly relativistic and an ISM mass of M /t] 
was collected. Since the reverse shock is just mildly relativistic at this stage, the shell's 
Lorentz factor have changed only by a factor of order unity. We call this the Newtonian 
Reverse Shock (NRS) case since the reverse shock is Newtonian relative to the unshocked 
shell. This is not to be confused with the fact that 7 ^> 1 in this case and the forward shock 
is ultra-relativistic. 

If initially £ < 1 then Rn is the smallest radius. The reverse shock becomes relativistic 
before it crosses the shell. At r > R N the reverse shock begins to reduce considerably 
the Lorentz factor of the shell's matter that it crosses. The Lorentz factor of the shocked 
material is 

7 (r) = /3/4 A -i/4 r -i/2 (6) 

(Sari, 1997). The shell has decelerated significantly by the time that the reverse shock has 
crossed the shell at Ra- At this stage, the Lorentz factor was reduced from its initial value 
of r] to r]£ 3 > 14 = (Z/A) 3/8 < 77. We call this the Relativistic Reverse Shock (RRS) case as 
most of the deceleration is done by a strong relativistic reverse shock. 
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The coasting radius, when the shell begins to coast freely, Rl, is related to the other 
radii by a simple relation, tj^ 2 Rl = (, 2 R S = V^Ra = -R 7 = Rn/C In the NRS case 
£ > 1 and Rl is the smallest radius. All the initial thermal energy is converted to kinetic 
energy before the shell begins to decelerate. In the RRS case (£ < 1) it is possible that 
the deceleration, due to external matter, begins before all the energy of the fireball was 
converted to kinetic energy. The conditions Rl < Rn and £ < 1 yield (//A) 3 / 8 < r] < \Jl/A 
(see figure [I|): £ should be larger than (A//) 1//6 in order that Rl < Rn- We limit the 
discussion to this case here,i?L < Rn, i.e., the fireball transforms all its internal energy to 
kinetic energy before the ISM begins influence its evolution. 

The NRS case is relevant if the initial fireball is small and contains relatively many 
baryons (equivalently, rj is relatively small), while the RRS case is relevant if the fireball 
is large and it is less polluted by baryon. The shell has given the ISM most of its energy 
either at R 1 {= Rn = Ra) in the NRS case, or at Ra in the RRS case, 

2.3. Relativistic Self Similar Deceleration 

For r > _R 7 for the NRS case, or r > Ra for the RRS case, the shocked ISM contains 
most of the energy E. From this stage the shell plays a negligible part in the consecutive 
evolution. The profile of the shocked ISM is determined now just by two parameters, E and 
Pi- 

When the forward shock reaches a radius R the fireball has a mass 47rpii? 3 /3. The 
energy in the shocked fluid is therefor oc p\R Zr y 2 . Since this equals the total energy of the 
system E, we obtain a scaling law, 7 oc (E / ' piY^R' 3 ^ 2 . The exact proportionality constant 
depends on the profiles behind the shock. Blandford and McKee (1976) describe an analytic 
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self similar solution, in which the Lorentz factor, the density and the pressure are given by 

7 (f , r) = ^fX~ l, \ P(t, r) = 2V2 Pl T X - 5/ \ p(t, r) = Ip^x' 17 ' 12 , (7) 

where T(t) = y / 17/87r(i?//)~ 3 / 2 is the Lorentz factor of the shock itself, and the similarity 
variable x is defined by x(t> r ) = 1 + 8r 2 (l — r/R). The shock radius is given by 
R = R(t) ~ (1 — l/8r 2 )t. In the previous section r denotes the radius of a fluid element in 
the shell, but r and t are independent coordinates here. 



2.4. The Sedov- Taylor Solution 

The Blandford-McKee self similar solution is derived with the assumption 7, T 3> 1. 
This assumption breaks down when the shock sweeps out a volume ~ Z 3 of ISM and its 
motion becomes non-relativistic. At the radius I the Sedov- Taylor (Sedov, 1946; Von 
Neumann, 1947; Taylor, 1950) solution becomes a good approximation . A characteristic 
length scale at a time t, which we can form from the two parameters E and pi, gives 
the shock radius R(t) = a(Et 2 / pi) 1 ^ within a numerical constant factor a depending 
on the adiabatic constant 7 (a — 0.99 for 7 = 4/3). The velocity of the shock is 
u = dR/dt = 2R/5t. The velocity /3 2 , the density p 2 an d the pressure p 2 just behind the 
shock can be expressed in terms of u: j3 2 = 6m/7, p 2 = 7pi,p 2 = 6piU 2 /7 for the adiabatic 
index 7 = 4/3. The profile throughout the region behind the shock is given by the velocity 
(3 = 2rV/5t, the density p = p±G and the pressure p = 3r 2 p 1 GZ/25t 2 (Landau & Lifshitz, 
1987). The dimensionless variables V, G, Z are functions of the similarity variable ( = r/R 
only and are given in an implicit analytical form by the following equations: 

/7 \ -2 f 7 1 -232/99 f7 .5/11 

C 5 = (~V) {^(5-3V)} {-(4V-3)} , (8) 

Z = 2 -V 2 (l-V)(4V-3y\ (9) 

I f 7 9/11 f 7 n 116/33 

G = ^-^{^V-S)} {^( 5 - 3 ^)} • ( 10 ) 
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The pressure ration p/p2 tends to a constant as ( — > 0, while /9//3 2 oc C> p/P2 oc C 9 in this 
limit. 

3. Numerical Results 

3.1. Initial Conditions 

We consider an initial uniform spherical fireball surrounded by a uniform cold ISM. 
The initial conditions are determined by four parameters: the total energy E, the baryonic 
mass M (= E/rj) and the radius R of the fireball and the ISM density p 1 . A "cold" ISM 
means that its pressure is negligible compared to the pressure behind the shocks throughout 
the whole evolutions. For convenience, we set the initial time as Ro rather than zero. We 
have chosen two sets of initial conditions to represent the the RRS and the NRS cases: 
E = 10 52 erg, p! = lproton cm" 3 , r] = 50 and R = 3 x 10 10 cm for the NRS case (f = 43), 
and E = 10 52 erg, pi = lproton cm -3 , rj = 10 4 and R = 4.3 x 10 9 cm for the RRS case 

(e = o.i). 

We assume that the fluid is described by a constant adiabatic index 7 = 4/3 although 
in reality this is not alway true. The shell's matter may cool during the coasting phase. 
However at the end of the coasting stage, the reverse shocks heat the fireball shell, and the 
ISM shocked by the relativistic forward shock carries most of the energy E. The main part 
of the system is subject to relativistic particles again. The forward shock is decelerated as 
T oc i?~ 3 / 2 and it becomes Newtonian at /. After that the scaling laws in the Newtonian 
regime depends on the adiabatic constant. But even then the shocked electrons will remain 
relativistic for a long time after this transition and the adiabatic index will remain around 
4/3. 



3.2. Numerical Results 



Figures depict the evolutions of the Lorentz factors for the NRS and the RRS cases, 
respectively. One clearly sees the initial acceleration phase in which the Lorentz factor 
increases linearly with time. This is followed by a coasting phase in which the Lorentz 
factor is a constant: rj. This stage ends when the effect of ISM becomes significant and 
most of the kinetic energy is dissipated at i? 7 or Then a self similar phase begins in 
which the Lorentz factor of the forward shock decreases like T oc R~ 3 / 2 . At I the solution 
becomes non relativistic and it turns into the Sedov- Taylor solution. A difference appears 
between the NRS and the RRS fireballs only in the energy transfer phase. As expected a 
sharp decrease in 7 is seen in the RRS case. A more gradual transition is seen in the NRS 
case. 

In the following sections we compare the numerical results with the analytic estimates. 
We examine the validity of the estimates of Rl, R s , R^ i? 7 , R^ and I as indicators of 
transition scales. 

3.2.1. Free Acceleration Stage 

A highly relativistic shell is formed after a short acceleration phase (see figure ^). The 
width of the shell is constant in this acceleration stage (figure |5|). The density and the 
pressure peak at the same position in the local fluid frame. In the observer frame, the outer 
part of the shell has a higher Lorentz factor and the density peak is wider than the pressure 
peak. The density peak (dotted line) moves ahead of the pressure peak (solid line) in the 
observer frame (figure |5|) . 

The average Lorentz factor of the each fluid element in the shell increase as the radius 
of the shell increases, 7 ~ r/A. However, the Lorentz factor of the outermost layers is 
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well above the average. This is results from the initial sharp edges of the fireball. The 
initial acceleration of the outermost layers depends on the steepness of the initial pressure 
distribution at the edge of the fireball. The initial step function distribution that we have 
chosen leads to large acceleration of the outermost layers. We can regard the thin region at 
the boundary of the fireball as expanding in the free expansion velocity, independent of the 
fireball thickness A. However, this fast layer is thin and its mass is negligible. Except of 
the evolution of the maximal Lorentz factor (at the outermost layers) the evolution of the 
bulk of the fireball is not effected by the choice of the initial steepness of the boundary. The 
initial conditions are washed out later, when the interaction with the ISM is significant, and 
even the maximal Lorentz factor is then independent of the initial conditions. Therefore we 
do not discuss this maximal value, but instead we consider the average Lorentz factor over 
the "uniform shell". We define the average value as: 

(/) = J fmr 2 dr/ J mr 2 dr, (11) 

where m = 7 {p + (3 + j3 2 )p} is the effective mass density in the observer frame and the 
integrals are defined from the origin to a radius at which the Lorentz factor takes the 
maximal value. The evolution of the average Lorentz factor (7) for the NRS and the RRS 
cases is shown in figure and ^ (thick line) respectively. The Lorentz factors increase 
linearly with time. The average internal energy and the mass density in the observer frame 
are shown in figure || Rl is a good indicator to the transition from the acceleration stage 
to the coasting stage. The mass density equals the internal energy density at ~ Rl (see 
figure |). 

3.2.2. Coasting Stage 

After the fireball uses up all the internal energy at Rl, it coasts with a Lorentz factor 
1] (see figures |] and |3[). In the RRS case the shell has a frozen radial profile through this 
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stage, while in the NRS case the frozen pulse approximation breaks down at R s and the 
shell begins to expand before the ISM has most of the system energy at R 1 (see figure ||). 
We can see the transition of the scaling laws of the density and the pressure at R s in figure 



3.2.3. Energy Transfer Stage: The RRS Case 

In the RRS case the reverse shock becomes relativistic at Rn before it crosses the 
shell at Ra- From this moment onwards the Lorentz factor of the coasting shell is reduced 
considerably after the passage of the reverse shock. Figure [j] shows the deceleration by a 
relativistic reverse shock. There are four regions in the figure: the ISM, the shocked ISM, 
the shocked shell and the unshocked shell, which are separated by the forward shock (FS), 
the contact discontinuity (CD) and the reverse shock (RS). 

Using the shocks' jump conditions and the equality of the pressure and the velocity at 
the contact discontinuity, we can estimate the Lorentz factor of the shocked region, Jcd, 
and the Lorentz factor of the reverse shock, Trs, for a planner geometry (Sari & Piran, 
1995), 

lCD = f 1/4 ^fu/2, T RS = lcD /V2, (12) 

where / = p±j p\. The quantities just ahead of (or just behind) the reverse shock are 
denoted by the subscript 4 (or 3). For spherical geometry, the pressure is not constant in 
the shocked region bounded by the two shocks. The ratio / in the above equations should 
be replaced with (p2/pz)f- However, in our simulations, P2/P3 is a factor of a few at most. 
Then, equation (|1"2D can roughly explain the numerical result. The time it takes for the 
reverse shock to cross a distance dx in the shell material is dt = dxj4y/J '/2. As the reverse 
shock compresses the shell material, the width dx becomes dx/2 after the shock passes 
through it. Then, the distance between the contact discontinuity and the reverse shock is 
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t/274"*/7 (oc t 2 ). Figure § shows the numerical result. 

The relativistic reverse shock passes through the shell as t 2 and it decelerates the 
coasting shell drastically. After this drastic deceleration, the shocked shell slows down as 
t- 1 ' 2 (fi gure 0) due to the pressure difference between p 2 and p%. At R& when the reverse 
shock crosses the shell, 7cd reaches effectively the value of 72 expected from the relativistic 
self similar solution and the profile of the shocked ISM region begins to approach the self 
similar one. The transition into the self similar solution is shown in figure 

At the beginning of the self similar deceleration stage, the density of the shocked shell 
is much larger than that of the shocked ISM. There is a large gap of the density at the 
contact discontinuity. However, the density perturbation, as we discuss in section |], does 
not effect 7 and p, and it does not propagate in the local fluid frame. Then, as the blast 
wave expands, it leaves the gap, and the ratio between p 2 and the density of the shocked 
shell damps. 

3.2.4- Energy Transfer Stage: The NRS Case 

In the NRS case (£ > 1), the shell begins to spread at R s and the value of £ decreases. 
At R 7 , a coincidence Ra = R*y = Rn happens. The shocked ISM becomes the main 
component of the system and the profile approaches the Blandford-McKee solution. At this 
stage the reverse shock is just mildly relativistic, the shell's Lorentz factor changes by a 
factor of order unity (see figure [TOD . 
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3.2.5. Relativistic Self Similar Deceleration Stage 

In the deceleration stage, the hydro dynamic profiles of the shocked ISM depends on 
only E and p\. The numerical values of the Lorentz factor 72, the density p 2 and the 



pressure P2 just behind the shock are compared with the self similar solution in figure |H 
The numerical result of Lorentz factor is consistent with the Blandford-McKee self similar 
solution within a few % difference in the relativistic regime (see figure |ll| (b)). Though 
the density and the pressure peaks are narrower than the velocity peak (see equation (0)) 
and the numerical errors are larger, the density and the pressure agree with the analytic 
estimates within ~ 20%. In the RRS case, 72 begins to satisfy the relativistic self similar 
scaling at R = 1.9i?A (circle 1 in figure [TT] (b)). The self similar solution and the equation 
(§) with a factor l/27r 1//4 which we neglected for simplicity, give the analytic estimate 
1.5-Ra- In the NRS case, 72 reaches within a 10% error line at R = 2AR 1 . 



Flow profiles are plotted as a function of the similarity variable \ in figure |12| In the 
top panel of the figure, the value of the Lorentz factor just behind the forward shock is 
~ 15, and the value drops to ~ 5 in the inner region. The numerical results agree well with 
the self similar solution. In the bottom panel of the figure, the value of the Lorentz factor 
just behind the forward shock is ~ 8, and the value drops to ~ 2 at the inner region where 
the numerical results deviate from from the self similar solution. The relativistic self similar 
solution is derived with an assumption that each fluid element is highly relativistic, but it 
is a good approximation still at 7 ~ 5. 



3.2.6. Transition to the Sedov-Taylor Solution 

The Lorentz factor of the forward shock decreases and it becomes non-relativistic 
around Z(~ 1.9 x 10 18 in this case). The scaling laws of the velocity j3 2 , the density p 2 and 
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the pressure p 2 also gradually shift from the Blandford-McKee solution to the Sedov- Taylor 
solution around I (see figure [11]) . The circles 2 (R c2 = 0.46/) and 3 (i? c3 = 0.61/) in 
figure |TT] (b) give a rough estimate of the radii where the relativistic self similar solution 
becomes invalid and the Newtonian self similar solution becomes valid (72 ~ 1.9 at circle 2, 
fa ~ 0.70 at circle 3). The relation between fa and R in the Newtonian self similar solution 
is already valid at circle 3, but the shock radius R is still proportional to time t at this 
stage and it is smaller than the radius expected by the Newtonian self similar solution. The 
relation between fa and t in the Newtonian self similar solution becomes valid at circle 4 
(R ci = 1.8/) where fa = 0.14. 

The transition of the profiles in the shocked ISM region, from the relativistic stage to 
the Newtonian stage, is shown in figure |13[ The jump condition for a strong shock gives a 
simple relation between the ISM density p\ and the shocked density p' 2 {= P2I2) measured 
in the observer frame, p' 2 = 72(472 + 3)px- If ISM have been swept up at a radius R, the 
thickness of the blast wave is approximately R/^2^12 + 3) oc R for the relativistic stage. 
The thickness is oc R for the Newtonian stage. As the blast wave expands, it becomes 



broader (figure |T3|a). The inner part of the distribution begins to approach the Newtonian 
self similar one Around R c3 (figure |T3"|b). The density profile approaches the one expected by 
the Newtonian self similar solution at R c % (see figure |T3"|c). The pressure profile approaches 
the Newtonian one more gradually and it still evolves after R c3 (figure p~3|d) . The velocity 
and pressure profiles are consistent with the analytic estimates within 10% level at R c ^. 



4. Evolution of Perturbations to Blandford-McKee Solution 

The simulations presented in the previous sections have shown that the hydrodynamical 
profiles of the shocked ISM approaches the relativistic self similar Blandford-McKee solution 
at the end of the energy transfer phase. It implies that the self similar solution is attractive 
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and stable. In this section we consider the evolution of spherical perturbations to the self 
similar profile. 

For an infinite uniform fluid there are two types of perturbations in the linear theory. 
One type is a sound wave propagating with the sound velocity relative to the fluid. The 
other is an entropy-vortex wave moving with the fluid. It is an entropy perturbation with 
no change of the pressure and it includes a spherical density perturbation. 

As an unperturbed initial configuration we take a self similar blast wave with a total 
energy E = 10 52 erg in a uniform ISM of the density 1 proton /cm 3 , at the moment where 
the Lorentz factor behind the shock is 7 = 4000 (t ~ 5.3 x 10 15 ). We then add a Gaussian 
perturbation Sq/q = 5 exp [— (x — Xo) 2 / ^x 2 }, with S — l,Xo — 3 and Ax = 0.5 to one of 
the hydrodynamical variables: the Lorentz factor, the density or the pressure, q stands 
for these variables and 5q for its perturbation. Figures 14* (a), (b) and (c) depict the flow 



profiles as functions of x f° r the perturbation of p, 7 and p, respectively. The dashed 
lines are the initial profile (self similar solution + perturbation), and the solid lines are at 
St = 4.3 x 10 14 sec later. 



The density perturbation (figure [14] (a)) does not effect other quantities (7 and p) and 
it does not propagate in the local fluid frame. We use the self similar solution to estimate 
the evolution of Xe of the fluid element in the perturbed region (or any other fluid element). 
Taking the derivative of Xe along its line of flow we get dxe/dt = 4(% e + l)/t so that a fluid 
element which had been at Xo a t a time t , will be at Xeif) ~ Xo(^Ao) 4 - The fluid element 
which had a density po at to, will have a density p e (t) = po(^Ao) _13//2 - The ratio between the 
density perturbation Sp e (t,r(t)) and the unperturbed value p e (t,r(t)) is constant in time, 
but the perturbation will depart from the forward shock as Xe oc t 4 and the perturbation 
becomes less important as p e / P2 oc t -5 . Our numerical results agree well with these scalings. 

A velocity perturbation (figure [13 (b)) induces pressure and density perturbations. 
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The pulse of the coupled perturbations is decomposed after a short time into two pulses: 
An outgoing compression pulse (outgoing motion corresponds to a decreasing x) an d an 
ingoing rarefaction pulse (a rarefaction wave and a shock wave). These waves propagate 
with the speed of sound in the local fluid frame 7^ = ^/4p e /3p e = 2 r y^ 2 x~ 1 ^ 12 /3. The speed 
of the outgoing and the ingoing sound wave in the observer frame are 7s = 4 72 3/ V 7/12 /3 
and 372 //2 x -5 / 12 /4, respectively. Assuming that the initial position xo of the perturbation at 
a time to is not too far from the forward shock, 7s ^> T, we get the position of the outgoing 
pulse x — Xo(t/to)~ A - The pulse reaches the forward shock at t — Xc/^o an d it boosts the 
forward shock. The position of the ingoing pulse cannot be expressed by a simple analytic 
formula, but it departs from the forward shock much faster than a density perturbation. 
After the out-going pulse boosts the forward shock, the flow profile is almost the self similar 
one. 

A pressure perturbation (figure [L4| (c)) also induces perturbations in 7 and p. The 
perturbations consist of three components, two propagating component and a standing 
density perturbation. After an outgoing compression pulse and an ingoing shock propagate, 
a density perturbation still stays at the position where the initial pressure perturbation was 
(in the local fluid frame). The ingoing shock leaves the forward shock quickly. The density 
perturbation departs from the forward shock as discussed earlier. The outgoing pulse is 
basically the same one appearing in the case of an initial velocity perturbation. The pulse 
boosts the forward shock and at this stage the profile is almost the self similar one. 

5. Conclusions 

We have explored numerically the evolution of a relativistic fireball interacting with a 
uniform ISM, through the stages of initial acceleration, coasting, energy transfer to the ISM 
and the deceleration. These calculations begin when the fireball is at rest. They follow the 
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acceleration to a relativistic velocity and the subsequent slowing down to a velocity far below 
the speed of light. These calculations span more than eight orders of magnitudes in the 
size of the fireball. The current analytic understanding of the fireball evolution can explain 
well our numerical results. Initially, the Lorentz factor increases linearly with the radius 
during the initial free acceleration stage. At R L the fireball has transfered all its initial 
radiation energy to kinetic energy and it coasts. Then the energy is transfered to the ISM. 
This takes place at i? 7 for the NRS case, or at for the RRS case. After that the shocked 
ISM carries most of the initial energy of the fireball. The profile of the shocked ISM is then 
described well by the relativistic self similar Blandford-McKee solution. The forward shock 
decelerates as T oc R~ 3 ^ 2 and it becomes Newtonian at I. After that the non-relativistic 
self similar Sedov- Taylor solution sets in. We have shown that the relativistic self similar 
solution is an excellent approximation, down to 7 ~ 5 and a reasonable approximation 
even down to 7 ~ 2. We have also examined the transition into the relativistic self similar 
solution, and the transition from this solution to the non-relativistic self similar solution. 

With the spherical code, we have shown that the hydrodynamical profiles of the 
shocked ISM approaches to the relativistic self similar solution at the end of the coasting 
stage, even though the initial conditions of the simulation do not have self similar profiles. 
It implies that the relativistic self similar solution is attractive and stable for a spherical 
perturbation. We have tracked the evolution of the spherical density, velocity and pressure 
perturbations to the relativistic self similar solution in order to verify the stability of this 
solution. 

The Blandford-McKee solution is the basis of much of the GRB afterglow theory. It is 
used to obtain an explicit expression for the radial profile during the deceleration stage. We 
have shown that this model is reasonable even down to 7 ~ 2. It is also valid even in the 
case of a radial inhomogeneity in the ISM. 
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Though the results presented in this paper are intended as a more general explanation of 
the strong explosion problem rather than a detailed fits of the GRB afterglow, comparisons 
of the numerically predicted light curves based on this computations in detailed realistic 
models with the recent and upcoming observations will enable us to determine free 
parameters in the GRBs model that cannot be calculated from first principles, such as the 
total energy of the system, the surrounding density, the fraction of the energy that is given 
by the shocks to electrons or to the magnetic field. 

S.K. gratefully acknowledges the support by the Golda Meir Postdoc fellowship. R.S. 
thanks The Clore Foundations for support. This work was supported in part by a US-Israel 
BSF grant 95-328 and by a NASA grant NAG5-3516. 



-22 - 



References 

Blandford,R.D. k McKee,C.F. 1976, Phys. of Fluids, 19, 1130. 

Daigne, F. k Mochkovitch, R. 1997, preprint. 

Goodman,J. 1986, ApJ, 308, L47. 

Kobayashi,S., Piran,T. k Sari,R. 1997, ApJ, 490, 92. 

Landau, L.D. k Lifshitz,E.M. 1987, Fluid Mechanics 2nd ed. (Pergamon Press), Chap. X. 
Katz, J., 1994, ApJ, 422, 248.. 

Meszaros,P., Laguna,P. k Rees,M.J. 1993, ApJ, 415, 181. 
Meszaros,P. k Rees,M.J. 1992, ApJ, 397, 570. 
Narayan,R., Paczyhski,B. k Piran,T. 1992, ApJ, 395, L83. 
Paczyhski,B. 1986, ApJ, 308, L51. 

Piran,T., Shemi,A. k Narayan,R. 1993, MNRAS, 263, 861. 

Rees,M.J. k Meszaros, P. 1994, ApJ, 430, L93. 

Sari,R. 1997, ApJ, 489, L37. 

Sari,R. k Piran,T. 1995, ApJ, 455, L143. 

Sari,R. k Piran,T. 1997, ApJ, 485, 270. 

Sedov,L.I. 1946, Prikl. Mat. i Mekh.,10, 241 

Shemi,A. k Piran,T. 1990, ApJ, 365, L55. 

Taylor G.I. 1950, Proc. Roy. Soc. London, A201, 159. 

Von Neumann J. 1947, Los Alamos Sci. Lab. Tech. Series, Vol 7. 



-23- 




Fig. 1. — Allowed parameter region for I = 1.9 x 10 18 cm (E = 10 52 ergs and p\ = 
lproton/cm 3 ). The line (£ = 1) separates the NRS case (lower left) and the RRS case 
(upper right). R L < R N is the lower left region of the line (Rl = Rn)- 




Fig. 2.— 7 vs Time, for the NRS case (^ = 43, E = 10 52 erg, r) = 50 R = 3 x 10 10 cm). The 
average value (thick solid line), the value just behind the forward shock (thin solid line), the 
maximal value (dotted line) and the analytic estimate (dashed dotted line). 
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Fig. 3.— 7 vs time for the RRS case (£ = 0.1, E = 10 52 erg, r] = 10 4 #0 = 4.3 x 10 9 cm). 
The average value (thick solid line), the value just behind the forward shock (thin solid line), 
the maximal value (dotted line) and the analytic estimate (dashed dotted line). 
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Fig. 4. — Profiles of the internal energy density 4py 2 (solid line) and the mass density 
fry 2 (dotted line) in the frame of ISM during the initial acceleration stage (internal energy 
dominated stage to matter dominated stage). The initial parameters are similar to those in 
figure |) 
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Fig. 5. — The shell structures at different times in the acceleration and coasting stage. The 
internal energy density 4py 2 (solid line) and the mass density fry 2 (dotted line) normalized 
at the peaks. The initial parameters are similar to those in figure |2|) The initial time is set 
as R (= 3 x 10 10 ). 
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Fig. 6. — The average internal energy density 4py 2 (thick solid line), the average mass 
density pj 2 (thick dashed line) and an analytical estimate of the internal energy density 
(thin solid line) and the mass density (thin dashed line). The latter are normalized at the 
cross point of the numerical results. The initial parameters are the same as in figure |2|. 
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Fig. 7. — : The Lorentz factor 7, density p and pressure p in the deceleration stage by the 
reverse shock (t — 3 x 10 15 ) for an RRS solution. The x axis is the distance from the contact 
discontinuity (dashed line). RS and FS indicate the positions of the reverse shock and the 
forward shock, respectively. The initial parameters are the same as in figure [| 
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Fig. 8. — The distance between the contact discontinuity and the reverse shock as a function 
of time: r C o — r RS (solid line), tjl^^^ff (dashed line) for a RRS solution. The initial 
parameters are the same as in figure |3|. 
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Fig. 9. — The transition into the Relativistic self similar solution in a RRS solution. Profiles 
of 7, p and p at different times (t = 0A1R A , 0.6C1R A , 0.73-R A , 0.85-R A , 0.93i? A , Ra, 1-1-Ra 
and 1.2i? A ) The x axis is the distance from the contact discontinuity (dashed line). 



-31 - 



10 2 




Fig. 10. — The transition into the Relativistic self similar solution for a NRS solution. 
Profiles of 7, p and p at different times (t = 0.70i? 7 , 0.77i? 7 , 0.82i? 7 , 0.87i? 7 , 0.92i? 7 , 
0.98i? 7 , Fly, and l.li? 7 ) The x axis is the distance from the contact discontinuity (dashed 
line) . 




Fig. 11. — The deceleration stage. (a)/3272 (solid lines), density (dashed lines) and pressure 
(dashed dotted lines) just behind the shock vs the shock radius. The numerical results are 
thick lines and the relativistic and the Newtonian self similar solution are thin lines. We plot 
72 (A) rather than /?272 for the relativistic self similar solution (the Newtonian self similar 
solution). The initial parameters are the same as in figure |3|. 

(b) The ratio between the numerical results and the analytic estimates. 72 or (3 2 (solid line), 
density (dashed line) and pressure (dashed dotted line) just behind shock and the shock 
radius (dotted line). The numerical results are compared with the relativistic (thick line) or 
the Newtonian (thin line) self similar solution. 72 (fa) is compared with the analytic estimate 
in the relativistic regime (the Newtonian regime). The circles 1,2,3 and 4 indicating the cross 
points are at R = 1.9R&, 0.46/, 0.61/ and 1.8/ respectively. 




X 

Fig. 12. — The Lorentz factor (solid lines), the density (dashed lines) and the pressure 
(dashed dotted lines) as a function of x during the relativistic deceleration stage at different 
times: (a)t = 2 x 10 17 (72 ~ 15) (b)t = 4 x 10 17 (72 ~ 8): The numerical results are thick 
lines and the analytical results are thin lines. 7, p and p are normalized with the numerical 
results 72, P2 and p2- The vertical lines indicate the positions where the fluid elements have 
the values of the Lorentz factor 7 = 2, 5 or 10. The initial parameters are the same as in 
figure §. 
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Fig. 13. — The transition of the flow profile of the shocked ISM from the relativistic 
Blandford-McKee solution to the Newtonian Sedov- Taylor solution: (a) (3 vs r (b) 6/62 
vs r/R (c) p/p2 vs r/R and (d) p/p2 vs r/R. The dotted lines shows the Newtonian self 
similar solution. The number denote the shock radius R =: 1) 0.21/, 2) 0.33/, 3) 0.52/, 4) 
0.76/, 5) 1.0/, 6) 1.2/, 7) 1.3/, 8) 1.5/, 9) 1.6/, 10) 1.8/, 11) 2.0/, 12) 2.2/, 13) 2.4/. The initial 
parameters are the same as in figure |3|. 
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Fig. 14. — Stability of the Relativistic self similar solution. Perturbation of p - (a), 7 - 
(b) and p - (c). Shown are the unperturbed self similar solution (dotted lines), the initial 
perturbed profiles (dashed lines) and the future evolution (solid lines) 



